Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 19 October 2009 (MN WTjiK style file v2.2) 



New Optimised Estimators for the Primordial Trispectrum 



Dipak Munshi 1 ' 2 , Alan Heavens 1 , Asantha Cooray 3 , Joseph Smidt 3 , Peter Coles 1 , Paolo Serra 3 

Scottish Universities Physics Alliance (SUPA), Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK 
o School of Physics and Astronomy, Cardiff University, CF24 3AA 

Department of Physics and Astronomy, University of California, Irvine, CA 92697 

CN ■ 



O 

u 

43 
6 



> 

CO 

m 



October 2009, Revision: 0.9 



ABSTRACT 

Cosmic microwave background studies of non-Gaussianity involving higher-order multispectra can dis- 
tinguish between early universe theories that predict nearly identical power spectra. However, the recov- 
ery of higher-order multispectra is difficult from realistic data due to their complex response to inho- 
mogeneous noise and partial sky coverage, which are often difficult to model analytically. A traditional 
alternative is to use one-point cumulants of various orders, which collapse the information present in a 
multispectrum to one number. The disadvantage of such a radical compression of the data is a loss of 
information as to the source of the statistical behaviour. A recent study by Munshi & Heavens (2009) 
has shown how to define the skew spectrum (the power spectra of a certain cubic field, related to the 
bispectrum) in an optimal way and how to estimate it from realistic data. The skew spectrum retains 
some of the information from the full configuration-dependence of the bispectrum, and can contain all 
the information on non-Gaussianity. In the present study, we extend the results of the skew spectrum 
to the case of two degenerate power-spectra related to the trispectrum. We also explore the relationship 
of these power-spectra and cumulant correlators previously used to study non-Gaussianity in projected 
galaxy surveys or weak lensing surveys. We construct nearly optimal estimators for quick tests and gener- 
alise them to estimators which can handle realistic data with all their complexity in a completely optimal 
manner. Possible generalisations for arbitrary order are also discussed. We show how these higher-order 
statistics and the related power spectra are related to the Taylor expansion coefficients of the potential in 
inflation models, and demonstrate how the trispectrum can constrain both the quadratic and cubic terms. 
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INTRODUCTION 

.he inflationa ry paradigm which solves the flatness, horizon and monopole problem makes clear predictions about the generation and nature of density 
'•perturbations iC-uthll 1 98 ll ; IStarobinskvll 19791 ; iLinddl 19821 : lAlbrecht & Steinhardtll 19821 ; ISatdl 198 lh . Inflationary models predict the statistical nature of 
these fluctuations, which are being tested against data from a range of recent cosmological observations, including the recently-launched all-sky cosmic 
microwave background (CMB) survey Planclfl Various ground-based and space-based observations have already confirmed the generic predictions of 
inflation, including a flat or nearly flat universe with nearly scale-invariant adiabatic perturbations at large angular scales. Several planned missions are 
also targeting the detection of the gravitational wave background through polarisation experiments - another generic prediction of inflationary models. 
The other major prediction of the inflationary scenarios is the nearly Gaussian nature of these perturbations. In the standard slow-roll paradigm, the 
scalar field responsible for inflation fluctuates wi th a minimal amount of self interaction which ensures that any non-Gaussianity generated during the 
inflat ion th rough self-interaction would be sma ll dSalopek & B on Jl 1 99Ct 1 1 99 ll ; iFalk etaill 19931 ; iGangui et al.ll 19941 ; lAcquaviva et al.|[2003l : iMaldacenal 
|2003|) . See|Bartolo, Matar rese & Riottol fe006l) for a recent review and more detailed discussion. Any detection of non-Gaussianity would therefore 
be a measure of self-interaction or non-linearities involved, which can come from various alternative scenarios such as c urvaton mechanism, warm 
inflation, ghost inflation as well as string theory inspired D-cc eleration and Dirac Born Infield (DBI) models of inflation ( Linde & Mukhanovl 1 19971 ; 

iGupta. Berera & Heaven sll2002l ; lLvth. Ungarelli & Wandsll2003h . 

Early observational work on detection of primordial non- Gaussianity, from COBE dKomatsu etai] |2002|) and MAXIMA JSantos et alj l2003h 
was followed by much more accurate analysis with WMAlQ Kom atsu et alj 1200 3l : ICreminelli et al.1 12007|; ISpergel et al] |2007|). Optimised 3-point 
estimators were introduced by iHeavenj jl998l). and have been suc c essively developed (jKomatsu, Spergel & Wandeltl 120051 : ICreminelli et al.1 120061 ; 
Creminelli, Senatore, & Zaldarriagdl2007l ; ISmith. Zahn & Dorell2007l : ISmith & ZaldarriagalbOOdT Indeed, now an estimator for Jnl which sat urates 
the Cramer-Rao bound has been found, capable of treating partial sky coverage and inhomogeneous noise dSmith. Senatore & Zaldarriaga 2009). The 



1 http://www.rssd.esa.int/index. php?project=Planck 

2 http://map.gsfc.nasa.gov/ 
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recent claim of a detection of non-Gaussianity in WMAP data ( Yadav & Wandelt 2008J) has given a tremendous boost to the study of primordial non- 
Gaussianity, as it can lift the degeneracy between various early-universe theories which predict nearly the same p rimordial power spectrum. Most de- 
tection strategies focus on lowest order in non-Gaussianity, i.e. the bispectrum or three-point correlation functions dKomatsu. Spergel & Wandelt|l2005l : 
ICreminellil2003l ; ICreminelli et al]2006l ; lMedeiros & Contaldcll2006l : ICabella et al. 12006 - lliguori et al]2007l ; ISmith. Senatore & Z aldarriaga l2009l) . This 
is primarily because of a decrease in signal-to-noise as we move up in the hierarchy of correlation functions - the higher-order correlation functions are 
more dominated by noise than their lower-order counterparts. Another related complication arises from the necessity to optimise such estimators, and the 
impact of inhomogeneous noise and partial sky coverage is always difficult to include in such estimates. 

The recent study by Muns hT& Heavensl d2009h suggested the possibility of finding optimised cumulant correlators associated with higher-order 
multi-spectra in the context of CMB studies. These corre lators are well-studied in the context of pr ojected surveys such as projected galaxy surveys or in 
the context of weak lensing studies using simulated maps dMunshil2000] ; lMunshi & Jairj 2000. 2001). Early studies involving cumulant correlators focused 
mainly on understanding gravity-induced clustering in collisionless media and were widely employed in many studies involving numerical simulations. 
Cumulant correlators are multipoint correlation functions, collapsed to two points Although they are two-point correlators, they carry information on the 
corresponding higher-order correlation functions. Due to their reduced dimensionality they do not carry all the information that is encoded in higher-order 
correlation function, but the y carry more than their one-point counterparts, namely the moments of the probability distribution function which are often 
used as clustering statistics dBernardeau et al.1120021) . 

One of the reasons to go beyond the lowest order in non-Gauss ianity was pointed out by m any authors, including, e.g., Riq uelmel & Sperge] 

(2006). At smaller angular scales, the secondary effect may dominate dSpergel & Goldberdll999al lbl), in direct contrast to larger angular scales where 
the anisotropies are generated mainly at the surface of the last scattering. These secondary perturbations are produced by interaction of CMB photons 
at much lower redshift with the intervening large-scale matter distribution. Such effects will be directly observable with Planck. The deflection of CMB 
photons by the large-scale mass distribution offers the possibility of studying the statistics of density perturbations in an unbiased way and provide clues 
to growth of structure formation for most part of the cosmic history. Weak lensing of the CMB can provide valuable information for constraining neutrino 
mass, d ark energy equation of stat e and also has the potential to assist detection of primordial gravitation waves through CMB polarisation information; 
see e.g. iLewis & Challinorl d2006h for a recent review. However weak lensing studies using the CMB need to address the contamination produced 
by other secondaries such as the thermal Sunyaev Zeldovich (tSZ) effects and kinetic Sunayev Zeldovich effect (kSZ) as well as by point sources , 
Although it is believed that these contaminations are not so important in case of polarisation studies. It was pointed out in iRiauelmel & Spergell J200^ > 
that a real space statistic such as {ST i (fl)ST(fl')) c (which is a cumulant correlator of order four) can be used to separate the kSZ effect or Ostriker- 
Vishniac (OV) effect from the lensing effect as the lensing contribution cancels out at the lowest order. It was shown that, in addition to quantifying 
and controlling the kSZ contamination of lensing statistics such statistics could as well play a very important role in providing new insight into the 
history of reionization. In a complete ly different context it was shown that this e stimator also has use for testing models of primordial non-gaussianity 
using redshifted 21cm observations dCooraviboOrj ; ICoorav. Li & Melc hiorri 2008). While cumulant correlators at third order {5T 2 (£l)8T(&')) c can 
provide information regarding the non-Gaussianity parameter /nl, fourth-order statistics such as (ST 2 (tl)ST 2 (tl')) c can go beyond the lowest order 
non-Gaussianity by putting constraints on the next-order parameter (?nl (to be introduced in later sections), albeit at lower signal-to-noise. However with 
ongoing CMB missions such as Planck the situation will improve and developing optimal methods for such higher-order cumulant correlators is the first 
step in this direction. There are now several studies which provide independent estimates of /nl however we still lack such co nstraints for anl- This 
clearl y is related to the fact that in typical models (Jnl is expected to be small, #nl < r / 50, where r is the scalar-to-tensor ratio dSerrv. Lidsev & Slothl 
2008). We also note that various studies have pointed out the link between the non-Gaussianity analysis and the estimators which test the anisotropy of the 
primordial universe. We plan to address these issues in a related publication. Several inflation ary models provide direct consistency relations between /nl 
and (? nl, e.g. gNL = (6/nl/5) 2 (in some publications it is also denoted as /2 or tnl e.g. (Hu & Okamoto 2002; Coorav 2006; [Serrv. Lidsev & Slotr] 
2008)). Testing of these consistency relations can give valuable clues to the mechanism behind the generation of initial perturbations. However, a difficulty 
for methods designed to detect non-Gaussianity in the CMB is that other proc esses can contribute, such as gravitational le n sing, unsubtracted point s ources, 
and imperfect subtraction of galactic foreground emission discussed by, e.g jGoldberg & Spergell dl999h ; ICooravl (2000); Ve rde & Spergell 0002); Castro 
d2004l) ; lBabich & Pierpaohl d2008h . 

While the main motivation in this work is to study the primordial trispectrum, we note that mode-mode coupling resulting from weak lensing of 
the CMB produces a trispectrum which has been studied using power-spectra associated with (6T(Q) 2 ST(tl') 2 ) c . Lensing studies involving the CMB 
can achieve higher signal-to-noise ratio at the level of the bispectrum if we use external data sets to act as a tracers of large-scale structure. However, the 
lowest order at which the internal detection of CMB lensing is possible is the trispectrum. There ha ve been some attem pts to detect non-gaussianity at the 
level of th e trispectrum usin g e.g. COBE 4yr data release JKunz et. alj|200ll) . BOOMERanG data dTroia et. alj|2003l) and more recently using WMAP 3 
year data (Sper gel et al]2007l) . 

The layout of this paper is as follows. The section §2 we introduce the concept of the higher-order cumulant correlators and how they are linked 
to corresponding correlation functions in real space. In §3 we discuss the harmonic transforms of the cumulant correlators and their relations to the 
corresponding multi-spectra. We also discuss estimators based on pseudo-Ci (PCL) estimators used for power spectrum estimation and generalise them 
to two-point cumulant correlators at higher order in this section. In §4 we briefly discuss the "local" models for initial perturbations and the resulting 
trispectrum. These models are then used in §4 to optimise the power spectra associated with the multi-spectra. This approach is nearly optimal, describing 
the mode-mode coupling using the "fraction of sky" approach familiar from other studies. In §5 we present an estimator which is nearly optimal and can 
take into account partial sky coverage as well as realistic inhomogeneous noise resulting from the scanning strategy. This estimator can work directly 
with any specific theoretic al model for primordial trispectra, based on concepts of matched filtering, and generalises the results obtained previously by 
iMunshi & Heavensl d2009h . We present results both for one-point estimators and also generalise them to two corresponding power spectra. §6 is devoted 
to adding relevant linear correction terms to these estimators in the absence of spherical symmetry. Next, in §7 we introduce inverse covariance weighting 
and design a trispectrum estimator which is optimal for arbitrary scanning strategy. The estimator used in this section will also be useful also in estimating 
secondary non-Gaussianity. For homogeneous noise and all-sky coverage the estimators are identical. Finally, §8 is devoted to concluding discussions 
and future prospects. 
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2 CORRELATION FUNCTIONS AND THE CUMULANT CORRELATORS 

The temperature fluctuations of the Cosmic Microwave Background (CMB) are typically assumed to be a realisation of statistically isotropic Gaussian 
random field. For Gaussian perturbations all the information needed to provide a complete statistical description is contained in the power spectrum of the 
distribution; for a non-Gaussian distribution higher-order correlation functions are also needed. With the assumption of weak non-Gaussianity only the 
first few correlation functions are needed to describe the departure from Gaussianity. We will denote the n-point correlation function by £n(&i, ■ ■ ■ , fXi)- 
The n-point correlation functions are decomposed into parts which are purely Gaussian in nature and those which sig nify departures from G aussianity. 
These are also known as connected and disconnected terms because of their representation by respective diagrams; see Bernardea u et alj d2002h for more 
details. At the level of the four-point correlation function the corresponding connected part, denoted by the subscript (. . .} c ), fi4 can be defined as: 

/x 4 A; . . . A) = {ST(Cli) . . . ST(Cl±)) c . (1) 

The connected component of the four-point function will be exactly zero for a purely Gaussian temperature field. The Gaussian contribution on the other 
hand can be written as a product of two two-point correlation functions. So the total four-point correlation function can be written as the sum of the 
connected and the disconnected part: 

£ 4 (^i . . . n 4 ) = & A A)& A, n 4 ) + & A A)& A> ft 4 ) + 6 A A))& A A) + ^ A • • • (2) 

As we will see, the Gaussian part will add to the scatter associated with any estimator for the four-point correlation function. At lower level £2(^1, ^2) = 
fi2(Cli, CI2) and £3(^1 . . . CI3) = . . . CI3) and hence there are no disconnected parts. The number of degrees of freedom associated with higher- 

order correlations increase exponentially with its order. This is mainly due to the increased number of configurations possible for which one can measure 
a higher-order correlation function. The cumulant correlators are defined by identifying all available vertices or points to just two points. There are two 
such cumulant correlators which can constructed at the four-point level. 

€si = <5 3 rA)<5TA)}c; £22 = <S 2 tA)£T 2 A)}c- (3) 

These two degenerate sets of cumulant correlators carry information at the level of four-point, but they are essentially two-point correla tion functions and 
can be studied in the harmonic domain by their associated power spectra. The first, £31 was studied in the context of 21cm surveys bv lCooravl ( l2006h and 
ICooray, Li"& Melchiorri ( 2008) and the second £22, was shown to have the power to separate the lensing contribution fro m the Kinetic Sunyaev Zel dovich 
effect (Riquelmel & Spergel 2006). These are natural generalisations of their third-order counterpart recentl y studied bv lMunshi & Heavens! d2009l). who 
introduced their optimised form for direct use on realistic data. These were later used bv lSmidt et alj J2009") to estimate /W, and lCalabrese et al. I J2009h 
to study lensing-secondary correlations from WMAP5 data. 

Although there are various advantages of working in real space, it is often easier to work in the harmonic domain. The main motivation to work 
in the harmonic domain is linked to the fact that inflationary models predict a well-defined peak structure for the power spectrum. These structures are 
well-known diagnostics for constraining cosmology at various levels. This is true for higher-order multi-spectra as they also involve the effect of the 
transfer functions. Note that the noise in CMB experiments is typically assumed to be Gaussian and will therefore contribute only to the disconnected 
terms. 



3 FOURIER TRANSFORMS OF CUMULANT CORRELATORS AND THEIR OPTIMUM ESTIMATORS 

The real-space correlation functions are clearly very important tools which can be used in surveys with patchy sky coverage. However recent CMB 
surveys scan the sky with near all-sky coverage. This makes a harmonic-space description more appropriate as various symmetries can be included in a 
more straightforward way. For Gaussian random fields with all-sky coverage the estimates of various statistics are loosely speaking uncorrelated. Even 
for a non-Gaussian field, they are reasonably uncorrelated at different angular scales or at different I as long as the non-Gaussianity is weak. We begin by 
introducing the harmonic transform of the observed temperature map ST(Cl) (f2 = (6, <j))) for all-sky coverage: 

a lm = J dft^jPy^ft) = J dil8T{h)Y{> m {Q). (4) 

Realistically however, we will only be observing the part of the sky which is not masked by Galactic foregrounds. The window function w(Q) which we 
will take as a completely general window can be used to define what is known as the pseudo harmonics which we designate as 2; m . 

~a lm ee J dClw(Cl)5T(Cl)Y; m (Cl). (5) 

We follow this description for rest of this paper. Any statistic X obtained from the masked sky will be denoted as X and the estimated all-sky version 
will be denoted X. The ensemble averages of the unbiased all-sky estimators which coincide with the theoretical models will be denoted only by the 
corresponding latin symbol X. 

(3 1) (2 2) 

3.1 Two-Point Estimators for the Trispectrum C\ ' , C, ' 

At the level of four point cumulant correlators we have two different estimators which can independently be used to study the trispectrum. These estimators 
that we discuss here are direct harmonics transforms of these two point correlators. We introduced the cumulant correlators in Eq.([3}. The two estimators 
we define in this section are related to {ST 2 (Cl)S 2 T(Cl')) and {8T 3 (Cl)ST(Cl')) through harmonic transforms. For the cubic function ST 3 (Cl), we denote 

the harmonic transform as and similarly a^l is the harmonic transform of the quadratic form 5T 2 (£l). We obtain the following relations for all-sky 

(2) 

coverage, and for the pseudo-harmonics a lr ^ with a mask. 
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a T2, = X] ai i"n a <2™2 y rfOV' iimi (f2)Y i2m2 (r2)Y i ^(A); (6) 
= ^ • • • ^ a hmi ai 2m2 wi 3 , n3 J dClY limi ({l)Yi 2m2 (Cl)Yi 3m3 (Q,)Y* m (Cl) = ^ ^imi'm'a^/- (7) 

l^mi ^3 m 3 (I'm') 



- ; (2) 



The corresponding results for and 5^ are as follows: 



a ;m = X/ '" X/ a hm 1 ai 2 m 2 ai 3 m 3 J dClY hmi (Cl)Y hm2 (Cl)Y l3 , n3 (€l)Y* m (£l); (8) 
5 im = X/ "" X] a 'i">i a '2"i2 a i3™3™i4'"4 J dClY hmi (Cl)Y hrn2 (Cl)Yi 3m3 (Cl)Y u , n4 (Cl)Y* lm (Cl) = ^2 Kimi'm'al?^- (9) 

i]_rai Z47714 * (I'm') 

The coupling matrix A^ m z'm' encodes information about the mode coupling which is introduced because of the masking of the sky faivon et alfeOOlh : 

Ki imihm2 [w] = y w(fl)Y"i imi (n)yj 2m2 (fl)dfl = 



Wl ™[ ^ J ( J ( mr m 2 m 3 >' (10) 

'3^13 



The matrices here denote the 3J symbols tedmondsi ri968). Using the harmonic transforms we define the following power spectra which can directly 
probe the trispectra. 

c (2,2) _ 1 (2). (2). J_^-(2)..(2) Qn 

??i m 

From the consideration of isotropy and homogeneity we can write the following relations: 

\ a im a im ) c ^ L l ll'°mm', \ a im a lm) c — (--; 0n>O mm >. (IZ) 

(2 2) 

The pseudo-power spectrum which is recovered from the masked harmonics is defined in an analogous way. Finally the resulting power spectra Cy ' 
can be expressed in terms of the trispectrum T/ 1 / 2 (I) by the following expression: 



.(2,2) _ rhUm / (2fx + l)(g a + l) / (2Z 3 + l)(2l 4 + l) ( h h I \ h h I \ 

- 2^ J M 2 Wy 4tt(2Z + 1) V 4ir(2i + 1) 1 M ]■ 



The trispectrum which is in troduced h ere is a four-point correlat ion function in harmonic space. The definition here ensures that it is invariant under 
various transformations; see iHul d2000h and lHu & Okamotol j2002h for detailed discussion: 

(a limi a l2m2 a l3m3 a limi ) c ^^- l ) MT hi^ L ) ( ^ M ) [ m 3 m 4 —M ) ' (14) 

LA/ V 7 V 7 

Partial sky coverage can be dealt with in an exactly similar manner. By using the harmonic transforms of the masked sky and expressing the masked 

J~i ( 2 2) f 2 2) 

harmonics in terms of all-sky harmonics we can relate the PCL versions of these power spectra C{ ' in terms of its all-sky components Cy ' . This 
involves the harmonic transform of the mask Wi m and its associated power spectrum Wi. The coupling matrix for the M n i is the same as that which is 
used for inverting the PCL power spectrum to recover the unbiased power spectrum. 

(2,2) _ sr^^f 1 i' i" Visr + i), ,, 



/ 4tt 



ZlZ2 I3L 

where 



X V 7*Wz"l / (^l + l) (^2 + 1) (2^3 + 1) (2/4 + 1) !l Z 2 !' Z 3 i 4 l' \ _ ,(2,2) 

x 2. r i,««(Oy M 2F + i) V 4<2«' + l) I J I )-Z. M »' C '' ■ (15) 



2 



In the literature the power spectrum c\ ' is also referred as the second spectrum. The other degenerate cumulant correlator of the same order that contains 
information about the trispectrum can be written as: 

C (3,D = 1 y R j f ( 3 ). (i) 1 c -(3,i) = 1 y Real ( a W. a W } . (17) 
l 2/ -)- 1 ^ L J ' 2/ -|- 1 ^ L J 
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The resulting power spectrum which probes the various components of the trispectrum with different weights can be expressed as follows: 

r (3,D_ V- T W 2rn / (2ii + l)(2k + iy l (2L + l)(2h+T) ( h h L\(L h I \ , , m 

C ' ~ t f? L lal [ } V 4tt(2L + 1) V 4tt(2; + 1) ^0 Oj{o J' 

With partial sky coverage we can proceed as before to connect the PCL version of the estimator C\ ' to its all-sky analogue. The resulting estimators 
combine the values of various components of the trispectrum with different weighting. The associated power- spectrum will therefore have a different 
dependence on parameters describing the primordial trispectrum. 



x V T l ^(D / (gi + lXgfa + l) / (2L + 1) (2*3 + 1) Mi la L W L * 3 /' \ _ Vm r (3,i) riQ . 
x T < 3 < (L) V 47 r(2i + l) V 4,(2I'+ 1) 1 J I ) ~ 2^ M "' C '' ' (19) 

The links to real space cumulant correlators are same as their third order counterpart: 

{5 2 T(Cl)5 2 T(ti')) c = — V(2Z + l)Pi(cos(fi ■ (l'))d 2 ' 2) '; <<5 3 T(f2)<5T(f2'))c = — V(2Z + l)P;(cos(0 ■ fV))^ 3 ' 1 '. (20) 

47T 47T ^— ' 

( I 

Hence we can conclude that it is possible to generalise these results to an arbitrary mask with arbitrary weighting functions. The deconvolved set of 
estimators at order p, q can be written as follows: 

Cr = M-'Cf/". (21) 

This is one of the important results of this paper. The mask used for various orders p, q is the same, but with the increasing order the number of degenerate 
power-spectra that can be constructed from a multi-spectrum increases drastically. As we have seen at the level of bispectrum we can keep one of the 
triangle sides fixed and sum over all contributions from all possible configurations of the triangle. Similarly for the trispectrum we can keep one of the 
sides of the rectangle fixed, or one of the diagonals fixed, and sum over all possible configurations. The possibilities increase as we move to higher-order 
in multispectra. Another related complication would be from the number of disconnected terms which need to be subtracted as they simply correspond 
to Gaussian contributions. At the 4-poin t level we need to sub tract the disconnected pieces from the dominant Gaussian component of temperature 
fluctuations including the noise jHi]|2000l : lHu & Okamotdl2002h : 



G \l\l( L ) = (~l) h+h y/(2h + l)(2i 3 + l)C h C h SL S llh S l2l3 + (2L + l)C h C h [(-l) l2+h+L 5 llh 6 hl4 + 6 h iJ hh ] , (22) 

where Ci is the power spectrum including noise. 

For the all-sky case and if we restrict only to modes with ordering l\ < I2 < I3 < I a, the non-zero component corresponds to terms with L = 
or h = I2 = I3 = I4. However for arbitrary sky coverage, which results in mode-mode coupling, no such general comments can be made. Estimators 
developed here in their all-sky form are known in the literature, and the results here show how to generalise them for partial sky coverage. However the 
main interest in data compression lies in using optimal weights for the compression, which is model-dependent as the weights depend on the specific 
model being probed. 

(3 1) f3 1) (2 2) (2 21 

If we introduce the Gaussian contributions to C\ ' as GJ ' , and to Cy as G\ ) , in the presence of partial sky coverage we can write: 
Gf» =Y,M lV GfV; G< a ' a > =£M„,G<?' a >. (23) 

V V 
For all-sky coverage one can obtain the corresponding results by replacing T^ 1 '/ 4 2 (L) by G'* (L) in respective equations: 



r (3.«_ >T r*,«vn / («i + l)(2h + l) / (2L + 1) (2*3 + 1) ( h h L L h I , ... 
G i - 2^ G <3.a L )y 4^(2L + l) V 4^(2/ + 1) olio I '~ 4 ' 



r (2,2) _ v-^ Ws,Um / (2il + l)(2i2 + l) / (2*3 + l)(2U + l) h h I \ (Is U I \ „„ 

°< - 4tt(2* + 1) V 4tt(2Z + 1) loOoHoOOT 

To subtract the disconnected Gaussian part, we use simulations. These are constructed from the same power spectrum that describes the non-Gaussian 
maps. Noise realisations which describe the actual map are also introduced in the Gaussian maps and exactly the same mask is used. This will ensure that 
the estimator remains unbiased. If we denote the final estimators after the subtraction of Gaussian disconnected parts by -D 3 ' 1 and Dp' 2 - 1 we can write 
for all-sky coverage: 

£(2,2) ^ £(2,2) _ £(2,2) . £(3,1) ^ £(3,1) _ £(3,1) _ (2fi) 

In the presence of a mask, deconvolved estimators are related to the masked estimators by the mixing matrix M u i : 

£(2.2) = M^{c[J- 2) - G< 2 ' 2) ); A' 3 ' 1 ' = E M w^' 1] ~ G f 1] )- ( 2? ) 

w w 
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Practical implementation of these estimators can provide a quick sanity check of a pipeline design for non-Gaussian estimators. While it is useful to 
keep in mind that these estimators are sub-optimal they are nevertheless unbiased and, depending on various choices of /nl and gNL, they can provide an 
analytical basis for computation of the scatter and cross-correlation among various estimators associated with different levels of the correlation hierarchy. 

(2 2^ (3 1) 

The dependence of D\ ' and D\ ' on /nl and gNL are different and hence can be used to provide independent constraints on both of these 
parameters without using third-order estimators and might be useful for providing cross-checks. 



3.2 One-Point Unoptimised Estimators for the Trispectrum cf 4 ^ 

The one-point cumulants at third order can be written in terms of the bispectra as: 

H 3 = (ST 3 (ft)) = 1-J ST 3 (ft)dft = h hhh B hl2h . (28) 

Similarly the one-point cumulant at fourth order can be written in terms as the trispectra as: 

M 4 = {5T 4 (ft)) c = i- J 5T 4 (ft)dft = -L h hh L h HU LTl^(L). (29) 

L I1I2I3I4 

We can also define the Sn parameters generally used in the literature as: 

S3=M3; 54=^4 — 3^2- (30) 

Throughout this discussion we have absorbed the experimental beam bi — exp{— 1(1 + 1) /2a 2 } in the harmonics a; m unless it is displayed explicitly. 
Here a b = FWHM/^8Zn2 for a gaussian beam. Alternatively it is also possible to define by respective powers of [12 to make these one-point estimators 
less sensitive to the normalisation of the power spectra. In that case we will have S3 = fJ-a/fJ^. These moments at fourth order are generalisations of the 
third-order moments, used as a basis for the construction of estimators for /nl by introducing the optimal weights A(r, ft) and B(r, ft). 



4 MODELS FOR PRIMORDIAL NON-GAUSSIANITY AND CONSTRUCTION OF OPTIMAL WEIGHTS 

The optimisation techniques that we introduce in next section follow the discussion in Munshi & Heavens! J2009l) . The optimisation procedure depends on 
construction of three-dimensional fields fro m the harmonic components of the temperature fields ai m with suitable weighting with respective functions 
which describes primordial non-Gaussianitvl Yada v. Komatsu,& W andelt 2007). These weights make the estimators act optimum and the matched filtering 
technique adopted ensures that the response to the observed non-Gaussianity is maximum when it matches with primordial non-gaussianity corresponding 
to the weights. 

In the linear regime the curvature perturbations which generate the fluctuations in the CMB sky are written as: 

aim = 4n(-i) 1 J ^<&(k)Af (ft)Y !m (fe). (31) 

We will need the following functions to construct the harmonic space trispectra as well as to g enerate weights for construction of optimal estimators (for 
a more complete description of predicting trispectrum from a given Inflationary prediction see ( IhuI|200(]| : Ihu & Okamotdl2002h ) : 

a l (r) = - k 2 dkA,(k)ji(kr); (3 t (r) = - / fc 2 dfcP <#> (fc)A i (fc)j i (fcr); w (r) = - / k 2 dk^i{k)ji{kr) (32) 
*■ Jo * Jo * Jo 

In the limit of low multipoles where the perturbations are mainly dominated by Sachs-Wolfe Effect the transfer functions A; (k) takes a rather simple 
form Ai(fc) = 1/3 ji (fer,) where r* = (770 — Tjdec) denotes the time elapsed between cosmic recombination and the present epoch. In general the transfer 
function needs to be computed numerically. The local Model for the primordial bispectrum and trispectrum can be constructed by going beyond linear 
theory in the expansion of the Additional parameters /ml and gNL are introduced which need to be estimated from observation. As discussed in 

the introduction, gNL can be linked to r the scalar to tensor ratio in a specific inflationary model and hence expected to be small. 

4>(x) = $ L (x) + f NL (4>i(x) - ($i(x))) + giVL<£-!(x) + Iinl (*|(x) - 3($|(x))) + . . . (33) 

We will only consider the local model of primordial non-Gaussianity in this paper and only adiabatic initial perturbations. More complicated cases 
of primordial non- Gaussiani ty will be dealt with elsewhere. In terms of inflationary potential V(4>) associated with a scalar potential <f> one can express 
these constants as (Hu 2000): 



5 1 d 2 hxV(<t>) 25 1 

/nl = — To - F< xj5 i # NL — 



6 8ttG dd> 2 * 54 8ttG 



d 2 \nV{<t>) \ _ ( d 3 lnV{(t>) \ ( d In V{<fi) \ 
d<j)' s J V d<j) J 



(34) 



There are two contribution to primordial non-Gaussianity. The first part is parametrised by /nl and the second contribution is proportional to a 
new parameter which appears at fourth order whic h we denote by (?nl ■ From th eoretical considerations in generic models of inflation one would expect 
jnl < r/5 with r b eing the scalar to tensor ratio JSerrv. Lidsev & Sloth 20081). 
Following jHi]|2 000) we can expand above expression in Fourier space to write: 
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*a(fc) = J ^$4k + k 1 )<E.2(k 1 )-(2^) 3 5 D (k) j glp^( kl ) (35) 

*s(*) = y ^$L(k 1 )$ ! (k 2 )<E- i *(ki). (36) 
The resulting trispectra associated with these perturbations can be expressed as: 

r*(*i,Aa,*s,*4) = (*(ki)*(ka)*(ks)«(k4)>c = J p^fo(k x + k 2 - K)fo(k 3 + k 4 + K)T s (ki,k 2 ,k 3 ,k4,K). (37) 
where the T(ki, k2, k3, k4, K) can be decomposed into two different constituents. 

Ti 2) (k 1 ,k 2 ,k 3 ,k4,K) =4/| ri P (K)P (k 1 )P 4 ,(k 3 ) (38) 
Ti 3) (ki,k 2 ,k 3 ,k4,K) = 5 ivi,{P0(K)P (ki)P*(k 3 ) + P (K)P (ki)P*(k3)} (39) 
The CMB trispectrum now can be written as 



T Wi( L )= 4:fN-Lh hl2L h hUL J r 1 dr 1 r 2 dr 2 F L (r 1 ,r2)ai 1 (r 1 )l3i 2 (r 1 )a h (r 2 )l3i i (r2) 

+gm.hi 1 i 3 z,hi s i i L / r 2 drf3i 2 (r)f3i 4 (r)[n h (r)(3 h (r) + fi h (r){3 h (r)]. (40) 



For detailed descriptions involving polarisation maps see jHu & Okamotoll2002l : lHi]|200(]| : lKomatsu & Spergelll200ll : lKogo etai1l2006l) . The CMB 
bispectrum which describes departures from Gaussianity at the lowest level can be similarly written as: 

B hhh = 2f NL h hhh J r 2 dr[a h (r)f3 l2 (r)p l3 (r) + a l2 (r)/3 l3 (r)/3 h (r) + ai 3 (r)/3 h (r)^ h (r)] . (41) 
We have defined the following form factor to simplify the display: 

(42) 



(2/r+l) (2/ 2 + l) (2/ 3 + l) f h h I 



I:', 



4tt y 

The extension to order beyond presented here to involve higher-order Taylor coefficients may not be practically useful as the detector noise and the 
cosmic variance may prohibit any reasonable signal-to-noise. 



5 NEAR OPTIMAL ESTIMATORS FOR TRISPECTRUM 

We will develop two estimators for extraction of power spectra associated with the trispectrum ( Hu 200Q|: IHu & Ok amoto 2002b in this section /C 3 ' 1 ^ and 
K. 2 ' 2 . These estimators are optimised version of similar estimators considered in the previous section. The derivation here parallels that of the previous 
section where unoptimised versions of these estimators were developed which uses the PCL type estimators. We will start by introducing four different 
fields which are constructed from the temperature fields. These fields are defined over the observed part of the sky and and are constructed using suitable 
weights to temperature harmonics. These w eights are functions of radial distance ai(r), Pi{r), fli(r) so the constructed fields are 3D fields. This method 
follows the same technique as introduced by Komatsu, Spergel & Wan deld J2005h . 



A(r,Cl) = ^A lm Y lm (£l); B(r,ti)=^2Bi m Yi m (Ci); M(r, fi) = ^ M lm Y lm (Cl); (43) 

lm Im I m 

m ji jr fJ>lm {AA\ 

A lm — -pq-dlmt J^lm = -pq-CLlm', ivllm = —p^-aim- (44) 

lw Ci Ui 

We have absorbed the beam smoothing into the corresponding harmonic coefficients. In addition to the weighting functions we will also need to define 
overlap integrals Pz,(ri,r 2 ) which act as a kernel for cross-correlating fields at two different radial distances. Note that the overlap integral depends on 
the quantum number L: 

F L ( ri ,r 2 ) = ^ y'fc 2 dfcPi, <s ,(A ; )j i (fcr 1 )ji,(fcr 2 ). (45) 

In the following subsections we will use specific forms for the trispectra to construct optimal estimator. Although the estimators developed here are 
specific to a given model clearly for any given model for the trispectra we can obtain similar construction in an optimal way. The weighting of these 
harmonics make the estimator a match-filtering one which ensures maximum response when the trispectra obtained from the data matches with that with 
theoretical construct. Inverse variance weighting ensures that the estimator remains optimal for a specific survey strategy. 
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5.1 Estimator for £ (4) 

The one-point estimators are simpler though they compress all available information in one number. One-point estimators can also be computed directly in 
real space and hence are simpler to compute when dealing with partial sky coverage in the presence of inhomogeneous noise. Using notations introduced 
above, we can write the one-point cumulant at fourth order by the following expression: 

K. {4) =Afl,J rj J r 2 2 dr 2 F{n,r 2 ) J dClA(r u £l)B(n, Cl)B(r 2 , Cl)A(r 2 , Cl) + 2g m J r 2 dr J dClA(r,Cl)B 2 {r,Cl)M{r,Cl) (46) 
which can be written in a compact form 

/C (4) =4/£ L [rj [ r 2 2 dr 2 F(r 1 ,r 2 )(A(r 1 ,{l)B(r 1 ,n)A(r 2 ,Q)B(r 2 ,Q)) + 2g m [ r 2 dr(A(r, ti)B 2 (r, Q)M(r, fi)>. (47) 



As expected there are two parts in the contribution. The first part depends on two radial directions through the overlap integral F(ri, r 2 ); the second is 
much simpler and just contains one line-of-sight integration. The amplitude of the first term depends on /f| L and the second term is proportional to ^nl- 
Typically they contribute equally to the resulting estimate. 

5.2 Estimator for K.\ 3 ' 1) 

Moving beyond the one-point cumulants we can construct the estimators of the two power spectra which we discussed before, C ; ' and C) ' . The 
corresponding optimised versions can be constructed by cross-correlating the fields A(ri, Cl)B(ri, Cl)B(r 2 , fi) with B(r 2 , f2). Clearly in the first case 
the harmonics depend on two radial distances n , r 2 for any given angular direction: 

A( ri )B(n)B(r 2 )\i m = J A{r 1 ,Cl)B(r 1 ,Cl)B(r 2 , A) Y lm (Cl) dti; A{r 2 )\ lm = J A(r 2 , Cl) Yi m {&) dCl. (48) 
Next we compute the cross-power spectra J^ B ' A {ri,T2) which also depend on both radial distances ri and r 2 : 

J? B ^ A {r U r 2 ) = ^^Y / F Uri,r2)Re a l[{A(r 1 )B(r 1 )B(r 2 )} lm {A(r 2 )}; m ]. (49) 

m 

The construction for the second term is very similar. We start by decomposing the real space product A(r, tl)B 2 (r, Q) and M(r, ti) in harmonic space. 
There is only one radial distance involved in both of these terms. 

A(r)B 2 (r)\i m = [ [A(r,Cl)B 2 (r,Cl)]Y lm (h) dti; M(r,Cl)\ lm = [ M(r)Y lm (Cl) dCl. (50) 



Finally the line-of-sight integral which involves two overlapping contribution through the weighting kernels for the first term and only one for the second 
gives us the required estimator: 

^ (3ll) = 4/ n 2 ,| rUnf r 2 2 dr 2 J l AB2 > A (r 1 ,r 2 ) + 2g nl [ r 2 dr£? B2 < M (r) (51) 

Next we show that the construction described above does reduces to an optimum estimator for the power spectrum associated with the trispectrum. 
The harmonics associated with the product field A(r\)B(r\)B(r 2 ) can be expressed in terms of the functions a(r) and /3(r): 

A(ri)B(ri)B(r 2 )\i m = a hmi a hm2 ai 3m3 ai 1 {ri)P h (ri)/3 h (r 2 )g^l^ 2M g^l sm . (52) 

LAI lm,limi 

The cross-power spectra J AB ,A (ri, r 2 ) can be simplified in terms of the following expression: 

Ji AB ' A (ri,r 2 ) = -—. - ^2(-l) M ^2 {F L {r 1 ,r 2 )a h (r 1 )ai 2 (r 1 )a h (r 2 )f3i(r 2 )} ( (53) 

LM m 

Where the Gaunt integral describing the integral involving three spherical harmonics is defined as follows: 



,,„ / (2Zi + 1)(2Z 2 + 1)(2Z 3 + 1) ( h h Is U h h k . - h 

y hhh - Y 47r ^ J \ mi m 2 m 3 J " 1 ' 

The second terms can be treated in an analogus way and takes the following form: 

c ab ifj r j = —^^(-^'^{oti^r^ai^r^ai^)^^^ (55) 

LM m 

Finally when combined these terms as in Eq.d51l>, we recover the following expression: 

= 27Tl£E ft^ft^W^ (56) 



I \ I2 L 
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The estimator )Cr depends linearly both on f^ L and <7nl- In principle, we can use the estimate of /nl from a bispectrum analysis as a prior, or we can 
use the estimators 5, \ /C, 3 ' 1 ^ and /Q 3 ' 1 -' to put joint constraints on /jyx and gML- Computational evaluation of either of the power spectra clearly will 
be more involved as a double integral corresponding to two radial directions needs to be evaluated. Given the low signal-to-noise associated with these 
power spectra, binning will be essential. 

5.3 Estimator for /Cf ,2) 

In an analogous way the other power spectra associated with the trispectra can be optimised by the following construction. We start by taking the harmonic 
transform of the product field A(r, Cl)B(r, fl) evaluated at the same line-of-sight distance r: 

A(r,Cl)B(r,Cl)\i m = [ A(r)B(r) Y lm (Cl) dtl; B(r,Q)M(r,Q)\ lm = / ' B{r)M{r)Y lm (tl) dtl, (57) 



and contract it with its counterpart at a different distance. The corresponding power spectrum (which is a function of these two line-of-sight distances n 
and r2) has a first term 

jf^fan) = — ^(n, ra)A(n, tyB{n, Cl)\i m A(r 2 , Cl)B(r 2 , ft)|,* ra (58) 

m L 

Similarly, the second part of the contribution can be constructed by cross-correlating the product of 3D fields A(Cl, r\)B(Cl, n) and B(Cl, r-2)M(Cl, r%) 
evaluated at different n and ri, with corresponding weighting function Fi(ri, ri): 

Cf B ' BM (r) = ^^^A(r,tl)B(rMi™B(r,ti)M(r,ti)\t m (59) 

m 

Finally, the estimator is constructed as: 

ICf- 2) = 4/^ | rUn J rldr 2 J* B ' AB {r^r 2 ) + 2g nl J r 2 dr£? B > BM {r) (60) 
To see they do correspond to an optimum estimator we use the harmonic expansions and follow the same procedure outlined before: 

Cf B ' BM {r) = ^J— ^(-l) m ^2 {a h (r)ai 2 (r)/3i 3 (r)^i i (r)} {ai imi ai 2m2 a hm3 a Umi )g™{™ 2m g™?J 4m . (61) 



j i AB > AB (r lj r 2 ) = -^—^{-l) M ^{F L {ri,r 2 )ai l {r 1 )pi 3 (ri)ai 2 ^^^ (62) 

LM m 

Here we notice that J AB - AB ( ri ; r2 ) ; s invariant under exchange of n and r2 but J AB - BM ( ri j r2 ) i s no t. Finally, joining the various contributions to 
construct the final estimator, as given in Eq.J60t. which involves a line-of-sight integration: 

£(2,2) = _j_ \p 1 t/ 1 ,' 2 (Or/ 1 / 2 (0- (63) 

The prefactors associated with f^ L and <?nl are different in the linear combinations K-'f' 2 ^ and /C» 3 ' , and hence even without using information from 

third-order we can estimate both from fourth order alone. 

Figs.Q]and[2]shows the primordial /C ; ' 2 ' 2 ' and /C; 3 ' 2 ' for the WMAP5 best-fit cosmological parameters dDunklev et al.l20 09). integrated over a range 
of 500 Mpc around recombination, and shown as a function of harmonic number I. The computations of ICi typically scales as l^ax w i m resolution which 
makes it quite expensive for high resolution studies. The computations are largely dominated by computations of 3J symbols. For a given configuration 
number of no n-zero 3J symbols for which the computations are required roughly scales as l^ax which explains the scaling. For more details about noise 
and beam see ISmidt et al.l J2009l) 



(2 2) I II L 

The unbiased version of /C, ' has been used in the context of study of lensing effects on CMB maps ( Coorav & Kesden 2003). While cross- 
correlational analysis can be helpful for detection of lensing on CMB maps for internal detection involving only CMB maps an analysis at the trispectrum 
level is necessary. 

5.4 Linking various Estimators 

(3 1) (2 2) 

The one-point estimators can be expressed as a sum over all Is of the estimators K\ ' ' or K.) ' . Therefore the corresponding optimised one-point 
estimators can be written as: 

f ^ _ T/ 3 / 4 (L)f/ 3 / 4 (L) 

£ (4) =E r r r r ' (64) 

Here 4 (L) denotes the direct estimator from the harmonic transforms and T^ 4 (L) act as weights which are constructed using the optimisation 
function discussed above. Similarly, at the two-point level the corresponding power spectra can be written as: 



10 Munshi et al. 



(21+ OK/ 2 '' 

4x10"' I : 




100 200 300 400 500 100 200 300 400 500 



1 1 

(2 2) 

Figure 1. The trispectrum-related power spectrum /C ; ' is plotted as function of angular scale I. The left panel is for a noise free ideal all-sky experiment. The right panel 
correspond to WMAP-V band noise. Various curves correspond to a specific parameter combination of (/jvi ,9nl) as indicated (see text for details). 



(21+l)K/ ; 




Figure 2. Same as[T]but for the trispectrum-related power spectrum /Cj ' plotted as function of angular scale I for ideal (left panel) and WMAP V-band noise characteristics. 



T/ 3 / 4 mf/ 1 / 2 (l) r , n ^ T} 1 } 2 (L)f/ 1 / 2 (L) 

^( 2 . 2 ) _ \ '1'2 > ' '3'4 > 1 . ]C ' \ '3' V ; '3' » ; (65) 

Ci 1 Ci 2 Ci 3 Ci i ' 1 Ci t Ci 2 Ci 3 Ci 



LL 



Each contribution from a specific configuration formed by various values of U and L is weighted by the corresponding Ci to make the estimator optimal. 
For computation of the variances the fields are however considered as Gaussian, which should be a good approximation as the fields are expected to be 
close to gaussian. Note that the Cis here take contributions from both the signal and the noise terms i.e. Ci = Cf + , where Cf is just the theoretical 



(2 2") f3 1) (2 1) 

Both K,] ' and K,\ ' involve both parameters /„; and g n i and can be used for a joint analysis, along with <S; ' to put constraints from realistic data. 



expectation for primordial perturbations. 
Both /C; 2 ' 2 ' and /C; 3 ' 1 ' involve both para: 
It is easy to check that the following relationship holds: 

£ (4) = ^(2Z + l)7vf ' 2) = ^(21 + l)ivf (66) 

l i 

These results therefore generalises the results obtained in (Munshi & Heaven J2 009) in the context of bispectral analysis. 



5.5 Subtracting the Gaussian or disconnected Contributions 

In estimating the trispectrum we need to subtract out the disconnected or Gaussian parts JhuI|200(]| : IHu & Okamotcll2002l) . To do this we follow the same 
procedure but replacing the simulated non-Gaussian maps with their Gaussian counterparts. Both maps should be constructed to have the same power 
spectrum, and identical noise. If the noise is Gaussian it will only contribute only to the disconnected part. For construction of an estimator which is 
noise-subtracted we need to follow the same procedure above by constructing Gaussian maps A G (r, Cl), B G (Cl, r) etc. 



A G (Cl,r) 



Y lm (Q); 



B G (Q,r) = ^B G n Y lm {Q); M G (n,r) =^M^ Y im (f2); 



aG _ a lm n G 



nG _ Plm G 



M? = Mi 

±y± lm 



Ci 



G 



(67) 
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We illustrate this using IC\ 2 ; the analysis is very similar for the other estimator. We start by replacing the quantities I AB ' AB (ri, r-z) and £f B ' AB (r 1 , r 2 ) 



by their Gaussian counterparts, V\ 



A G B G ,A G B G 



(n, ri) and TZ 



A G B G ,B G M G 



(r): 



p A G B G ,A G B G (ri)ra) = J- J2 F L (r 1 ,r 2 )A G (r 1 ,£l)B G (r 1 ,Cl)\ lm A G (r 2 , Cl)B G (r 2 , Cl)\t 



21 + 



(68) 



These quantities are then used for the construction of an estimator, which in practice aims to estimate the Gaussian contributions to the total trispectra. 

nf BG ' BGM \r) = ^^A G (r,(l)B G (rM^B G (r,n)M G (r,Cl)\t n 



21 + 



I 



(69) 



m 1 1 m 1 

As before we combine these quantities to form the weighted Gaussian part of the contribution to the trispectrum: 



Or=±f nl JrUr, frUr 2 V(n,r 2 ) + 2 9NL J r 2 drll(r) = 
Subtracting the Gaussian contribution from the total estimator we can write: 



(I). 



(70) 



(71) 



A very similar calculation for the other fourth-order estimator Q l ' provides an identical result. After subtraction of the Gaussian contribution we can 
write: 



(3,1) _ r (3,l) n (3,l) 



(72) 



The corresponding one-point collapsed estimator has the following form, and the relationship between the various estimators remains unchanged. 



^ (4) -EE 



-Tl£{l){Tlllt[L)-G\^{L)}. 



(73) 



In the next section we consider partial sky coverage and the resulting corrections to the Gaussian and non-Gaussian estimators. 



6 PARTIAL SKY COVERAGE AND OPTIMISED ESTIMATORS 

The terms that are needed for correcting the effect o f finite sky coverage and inhomogeneous noise are listed below. These corrections are incorporated by 
using Monte-carlo simulations of noise realisations 1 Creminelli 2003; Babich & Zaldarriaga 2004l: lBabich. Creminelli & Zaldarriagal20 04: Babich M2005l : 
ICreminelli et al.l2 006. 2007; Creminelli, Se natore, & Z aldarriaga 2007). Whereas in case of bispectral analysis it is just the linear terms which are needed 
for corrections, for trispectral analysis there are both linear and quadratic terms. 

We will treat the general case in later subsections, but first we present results for the simpler case of homogeneous noise and for high wavenumbers 
where the density of uncorrelated states is modified by inclusion of the fraction of sky observed, f s ky 



6.1 Corrective terms for the Estimator fCf' 1 in the absence of spherical symmetry 

J~A\B\B 2 ,A 2 1 \ ffA\B\B 2 ^A 2 .7-Lm .7-Quadl 

i - 7 — |yi ~ 1 ' J 

J sky 



j-(AiB 1 )B 2 ,A 2 _j_ J-\^1H2/*>1,A2 _|_ jJ-V°l-°2/^U,"fi2 _|_ rj H"! ("2 ,"-11 _|_ J-"l"2\o\,A2 1 _|_ jrJi 
(A 1 B 1 B 2 ),A 2 r r A 1 (B l B 2 ,A 2 ) ^■B 1 (A 1 B 2 ,A 2 ) ^■B 2 (A 1 B 1 ,A 2 ) 



(A 1 B 2 )B 1 ,A 2 



r{B 1 B 2 )A 1 ,A 2 



T-Lin 1 

( _ 1 

J sky 

jQuad ' 

( _ 7 

J sky 

The expressions are similar for the other terms that depend on only one radial distance 



A 1 B 1 (B 2 ,A 2 ) n -A 1 B 2 {B 1 ,A 2 ) ^B 1 B 2 (A 1 ,A 2 ) 



(74) 
(75) 
(76) 
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£ABB,M 



T-Lin J 
^-Quad 



fsky 



-A(BB,M) 
-I 



+ 2£ 



B(A,M) 
I 



+ 2C 



B(AB),M 
I 



fsky 



£A(BB),M _|_ 2£ B {AB),M 



(77) 
(78) 
(79) 



To simplify the presentation we have used the symbol A(n, fl) = Ai; A(r, Q,) = A and so on. Essentially we can see that there are terms which are 
linear in the input harmonics and terms which are quadratic in the input harmonics. The terms which are linear are also proportional to the bispectrum of 
the remaining 3D fields which are being averaged. On the other hand the prefactors for quadratic terms are 3D correlation functions of the remaining two 
fields. Finally putting all of these expressions we can write: 



Kf' 1} = 4/nl 



JrUnJ 



r 2 dr 2 J, 



AB ,A 



I 



(ri,r 2 )+ I r z drCf B M (r). 



(80) 



From a computational point of view clearly the overlap integral F^{r\,r 2 ) will be expensive and may determine to what resolution ultimately these 
direct techniques can be implemented. Use of these techniques directly involving Monte-Carlo numerical simulations will be dealt with in a separate 
paper (Smidt et al. in preparation). To what extent the linear and quadratic terms are important in each of these contributions can only be decided by 
testing against simulation. 

6.2 Corrective terms for the Estimator AT 2 ' 2 in the absence of Spherical symmetry 

The unbiased estimator for the other estimator can be constructed in a similar manner. As before there are terms which are quadratic in input harmonics 
with a prefactor proportional to terms involving cross-correlation or variance of various combinations of 3D fields and there will be linear terms (linear 
in input harmonics) with a prefactor proportional to bispectrum associated with various 3D fields. 



jLin 



A 1 B 1 ,A 2 B 2 _ 1 

fsky 

1 

fsky 



J~A 1 B 1 ,A 2 B 2 r Lin rQuad 
I -1, -1, 



jA 1 B 1 ,(A 2 B 2 ) + j 



(A 1 B 1 ),A 2 B ] 



_|_ jA 1 {B 1 ,B 2 )A 2 + jB 2 (A 2 ,B 2 )A 2 + jB 2 (A 2 ,A 2 )B 2 + jA 2 (B 2 ,A 2 )B 2 



^-Quad 1 

fsky 

The terms such as K J , 



J, 



A 1 {B 1 ,A 2 B 2 ) 



+ J,. 



B 1 {A 1 ,A 2 B 2 ) 



+ Jl 



(A 1 B 1 ,A 2 )B 2 



+ Jl 



(A 2 B 2 ,B 2 )A 2 



(81) 
(82) 
(83) 



(n , r2) can be constructed in a very similar way. We display the term AC, ' (n , r 2 ) with all its corrections included. 



fsky 



ZAB,BM r Lin r Quad 
— 1] — 1 i 



rQuad 



Jijin 



1 

fsky 
1 



C MB^),A + 2C , 



\D(B,A) _|_ j^B 2 {A,A) _|_ 2J£B{AB,A) ^{B 2 ,A> | 



fsky 



jA 1 (B 1 ,B 2 M 2 ) _,_ jB 1 (A ll B 2 M 2 ) + j 



(A 1 B 1 ,M 2 )B 2 _|_ j{A 2 B 2 ,B 2 )M 2 



(84) 
(85) 
(86) 



The importance of the linear terms greatly depends on the target model being considered. For example, while linear terms for bispectral analysis can 
greatly reduce the amount of scatter in the estimator for local non-Gaussianity, the linear term is less important in modelling the equilateral model. In 
any case the use of such Monte-Carlo (MC) maps is known to reduce the scatter and can greatly simplify the estimation of non-Gaussianity. This can be 
useful, as fully optimal analysis with inverse variance weighting, which treats mode-mode coupling completely, may only be possible on low-resolution 
maps. 



6.3 Corrective terms for the Estimator AC (4) in the absence of spherical symmetry 

The corrections to the optimised one-point fourth order cumulant can also be analysed in a very similar manner. It is expected that higher-order cumulants 
will be more affected by partial sky coverage and loss of spherical symmetry because of inhomogeneous noise. 



£< 4 > = 7 i-[/C 

J sky 



(A) 



(87) 



The linear and quadratic terms can be expressed in terms of MC averages. We list terms which involve the overlap integral and the one with single 
line-of-sight integration. 



4/nl / rfdri I r 2 dr 2 V F l L 2 {A 1 {B 1 A 2 B 2 ) + Bi(^i^ 2 B 2 ) + A 2 (AiBi^ 2 ) + B 2 {A 1 B 1 A 2 )} 



L 
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Table 1. Various multispectra and associated power-spectra are tabulated along with their cosmological use 



Order 


Cumulants & Correlator 


Power Spectra 


Use 


2 = (1 + 1) 


(5 2 T(n)), (ST(n)8T(n')) 


Ci 


Constraints on Cosmology (Q, Ho. . . .) 


3= (2 + 1) 


(8 s T(Cl)), (8 2 T(n)8T(n')) 




Inflationary Models /nl > Secondaries x Lensing bi 


4= (3 + 1), (2 + 2) 


(5 4 T(Q)), {8 :i T(Ci)S{n')),{S 2 T(n)8 2 T(n')) 


KW, k^ 2 \k^ 


/nl ■ 9nl i KSZ-Lensing Sep. Internal Lensing det. 



+2g NL J r 2 dr {A{B 2 M) + 2B(ABM) + M{AB 2 }} (88) 
7 Quad = 4/NL j r 2 dri J r2dr2 ^ F i2 {A 1 B 1 {A 2 B 2 ) + A 2 B 2 {A 1 B 1 ) + A 1 B 2 {A 2 B 1 ) + A 2 B 1 {A 1 B 2 )}) 

+2g NL J r 2 dr {AM{B' 2 ) + B 2 {AM) + 2BM{AB) } . (89) 
We have introduced the shorthand notation Fl 2 = i*i,(ri, ra)- 

7 REALISTIC SURVEY STRATEGIES : EXACT ANALYSIS 

In this section, we include the full optimisa tion, including the mode-mode coupling introduced by partial sky coverage and inhomogeneous noise , 
generalising results from the three-point level (Babich 2005; Smit h & Zaldamagal200d : ISmifh. Zahn & Dorel2007l : ISmith. Senatore &~Z aldarriaga 2009h . 
As before, we find that for the trispectrum, the addition of quadratic terms in addition to linera terms is needed. The analysis presented here is completely 
generic and will not depend on details of factorizability properties of the trispectra. For any specific form of the trispectrum the technique presented 
here can always provide optimal estimators. Importantly, this makes it suitable also for the study of the trispectrum contribution from secondaries, and 
offers the possibility of dete rmining whe t her an y observed connected trispectrum is primordial or not. Generalisation to multiple sources of trispectrum 
is straightforward, following lSmidt et al.l J2009h 

Trispectra from sec ondary anisotropies such as gravitational lensing are expected to dominate the contribution from primary anisotropies 
dCoorav & Kesdenll2003l) . The estimator we develop here will be directly applicable to data from various surveys, but the required direct inversion 
of the covariance matrix in the harmonic domain may not be computationally feasible in near future. Nevertheless an exact analysis may still be beneficial 
for low-resolution degraded maps where the primary anisotropy dominates. At higher resolution the exact analysis will reduce to the one discussed in 
previous sections. In addition it may be possible at least to certain resolution to bypass the exact inverse variance weighting by introducing a conjugate 
gradient techniques. 

The general theory for op timal estimation from data was developed bv lBabichl l l20olh for the analysis of the bispectrum, and was later implemented 
by Smith & Zaldarriaga (2006). For arbitrary sky coverage and inhomogeneous noise the estimator will be fourth order in input harmonics, and involves 
matched filtering to maximise the response of the estimator when the estimated trispectrum matches theoretical expectation. We present results for both 

(3 1 ) ( 2 2) 

one-point cumulants and two-point cumulant correlators or power spectra associated with trispectra. The estimators presented here, E\ ' and E\ ' are 
generalisations of estimators Cj 3,1 \ C^ 2 ' 2 ^ and IC^ 3 ' 1 ^ and /C' 2 ' 2 ' presented in previous sections. 



7.1 One-point Estimators 

We will use inverse variance weighting harmonics recovered from the sky. The covariance matrix, expressed in the harmonics domain, [m when used 
to filter out modes recovered directly from sky are a; m are denoted as Alm- We use these harmonics to construct optimal estimators. For all-sky coverage 
and homogeneous noise, we can we recover Ai m = ai m /Ci, with C;s including signal and noise. We start by keeping in mind that the trispectrum can 
be expressed in terms of the harmonic transforms a; m as follows: 

W^ + DEE^fi t £)(£ t - M ) a ^- a ^ (iZa, b ,c,d). (90) 

Based on this expression we can devise a one-point estimator. In our following discussion, the relevant harmonics can be based on partial sky coverage. 

LAI l imi x 7 v 7 

The term A(U, L) is introduced here to avoid contributions from Gaussian or disconnected contributions. A(U, L) vanishes if amy pair of ks becomes 
equal or L — which effectively reduces the trispectra to a product of two power spectra (i.e. disconnected Gaussian pieces). 

We will also need the first-order and second-order derivative with respect to the input harmonics. The linear terms are proportional to the first 
derivatives and the quadratic terms are proportional to second derivatives of the function Q[a], which is quartic in input harmonics. 
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LM lirrii V 7 V 7 

The first-order derivative term di m Q^ [a] is cubic in input maps and the second order derivative is quadratic in input maps (harmonics). However unlike 
the estimator itself [a], which is simply a number, these objects represent maps constructed from harmonics of the observed maps. 

^ m ^) H = I^(-l)-^A fe L)T^(L)(^ £ MP £ _ L M )a lama a hmb . (93) 

LM liTTli v 7 v 7 

The optimal estimator for the one-point cumulant can now be written as follows. This is optimal in the presence of partial sky coverage and most general 
inhomogeneous noise: 

E (4) H = i [Q^yC^a] ~ [C^aUid^Q^lA]) - [C- 1 a] tm [CT 1 a] Vml (ai m a l , m ,Q w [C- 1 a])} . (94) 

The terms which are subtracted out are linear and quadratic in input harmonics. The linear term is similar to the one which is used for bispectrum 
estimation, whereas the quadratic terms correspond to the disconnected contributions and will vanish identically as we have designed our estimators in 
such a way that it will not take any contribution from disconnected Gaussian terms. The Fisher matrix reduces to a number which we have used for 
normalisation. 

F = i= iEE E E (-i) M (-i) M 'A(/ i ;L)A(4iO^;/;(^;i , w 

LM L'M' (all Im) (all I'm') 

I l a h k l 'd L | r^ 1 r^ 1 (Q^\ 

\ m' a m' b M' )\m' c m' d —M' J ^'^a,i' a m' a ■ • • ^i d m d ,i' d m' d - V J > 

The ensem ble average of this one-point estimator will be a linear c ombination of parameters f^ L and #nl • Estimators constructed at the level of three-point 
cumulants dSmith & Zaldarr iaga 2006: Munshi & Heavens! 20091) can be used jointly with this estimator to put independent constraints separately on /nl 
and £nl • As discussed before, while one-point estimator has the advantage of higher signal-to-noise, such estimators are not immune to contributions 
from an unknown component which may not have cosmological origin, such as inadequate foreground separation. The study of these power spectra 
associated with bispectra or trispectra can be useful in this direction. Note that these direct estimators are computationally expensive due to the inversion 
and multiplication of large matrices, but can be implemented in low-resolution studies where primordial signals may be less contaminated by foreground 
contributions or secondaries. 



la lb L 

m a nib M 



lc Id L 

m c md —M 



7.2 Two-point Estimators: Power Spectra associated with Trispectra 

Generalising the above expressions for the case of the power spectrum associated with trispectra, we recover the two power spectra we have discussed 
in previous sections. The information content in these power spectra are optimal, and when summed for over L we can recover the results of one-point 
estimators. 

Ql 2) M = 77 E(-1) M T TlfAL) ( L k L M )( lc ld L M ) a lama . . . a ldmd . (96) 

M ijTOj v 7 v 7 

The derivatives at first order and second order are as series of maps (for each L) constructed from the harmonics of the observed sky. These are used in 
the construction of linear and quadratic terms. 

OimQ^H ^B-D-E^^W ( I ^ L M ) ( X ^ _ L M ) a lama ...a lemc . (97) 

T li-rrii V 7 V 7 

We can construct the other estimator in a similar manner. To start with we define the function Q^ 3,1 ' [a] and construct its first and second derivatives. 
These are eventually used for construction of the estimator E ( L 3 [a]. As we have seen, both of these estimators can be collapsed to a one-point estimator 
[a]. As before, the variable a here denotes input harmonics a; m recovered from the noisy observed sky. 

^H^EEM^-E^^^S^fM X s T )(t t 4) — --— <*> 

M ST lirrii V 7 V 7 

The derivative term will have two contributing terms corresponding to the derivative w.r.t. the free index {LAI} and the terms where indices are summed 
over e.g. {Im}, which is very similar to the results for the bispectrum analysis with the estimator Q^ 2 ' 1 ' [a]. One major difference that needs to be taken 
into account is the subtraction of the Gaussian contribution. The function A(h, L) takes into account of this subtraction. 

lm Qf [a]= EE T KW CD ( I tb t){ t td —T ) • • • W 

ST lirrii V 7 X 7 
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T 



+ 3 EE ffl «'E( M I Til 1 m, —T K»»<W ( 100 > 



L I S \ l b l c S 
M m T J \m b m c -T 

(3 1) (2 2) 

Using these derivatives we can construct the estimators E L ' and E L ' : 

1] = NZh [Q^yC^a] - [C7- 1 a] Jm <ft m Qg' i:, [Cr 1 a])} (101) 



El 



(2.2) 



Nil [Q^ 2) [C- l a] - [C- 1 o] im (a m Q^' 2) [C- 1 a])} , (102) 

where summation over 1/ is implied. The quadratic terms will vanish, as they contribute only to the disconnected part. The normalisation constants are 
the Fisher matrix elements Fll 1 which can be expressed in terms of the target trispectra T l L ^{L) and inverse covariance matrices C" 1 used for the 

(3 1) (3 1) 

construction of these estimators. The Fisher matrix for the estimator E L ' , i.e. Fi£> can be expressed as: 



[n-% l , = = (i) y y (-in-i) M xtw<r/>')] 

A(W;L)A(li;L') x 



ST,S'T' (all lm,l'm') 

la lb S \ ( L Id S \ ( l a 1^ S \ ( L 1^ S 



m a m b T J \ M m d -T J \ m' a m' b T' J \ M' m' d —T' 

X LM ,V M 1 ^ l a m a X^m'S" l b m b ,l' b m' b ^l t m c , l' c m' c + ^^LM ,l' a m' a ^l a ™a , L' M' ^ l b m b ,l' b m' b m c ,l' c m' c } ' (103) 

(2 2) (2 2) 

Similarly for the other estimator E L ' , the Fisher matrix F LL , can be written as a function of the associated trispectrum and the covariance ma- 
trix of various modes. Fo r further simplification of these expressions we need to make simplifying assumptions for a specific type of trispectra, see 
iMunshi & Heavens! J2009h for more details for such simplifications in the bispectrum. 



la lb L \ I L Id L 

m a m b M ) \ M m d -M 



[N-% L ,=F^ = Q YY, E E {-l) M i-l) M 'Tl±{L)Tlf«{L') 

M M' (all Im) (all I'm') 

x (i 4 m')(m' l X i#0 A(W;i)A( ^'^^^ (104) 

Knowledge of sky coverage and the noise characteristics resulting from a specific scanning strategy etc. is needed for modelling of C, _1 ,, , . We 

<d m d,' d m d 

will discuss the impact of inaccurate modelling of the covariance matrix in the next section. The direct summation we have used for the construction 
of the Fisher matrix may not be feasible except for low resolution studies. However a hybrid method may be employed to combine the estimates from 
low resolution maps using exact method with estimates from higher resolution using other faster but optimal techniques described in previous section. In 
certain situations when the data is noise-dominated further approximation can be made to simplify the implementation. A more detailed discussion will 
be presented elsewhere. 



7.3 Approximation to exact C 1 weighting and non-optimal weighting 

If the covariance matrix is not accurately known due to the lack of exact beam or noise characteristics, or due to limitations on computer resources, it can 
be approximated. An approximation 7? of C _1 then acts as a regularisation method. The corresponding generic estimator can then be expressed as: 

Ella] = E^ 1 ]"-' {QlW - [Ra]i m {di m Q z L ,[Ra]}} ; Z g {(2,2), (3, 1)}. (105) 

L' 

As before we have assumed sums over repeated indices and (■) denotes Monte-Carlo (MC) averages. As evident from the notations, the estimator above 
can be of type E^' 1 ^ or E^' 2 K For the collapsed case E^ can also be handled in a very similar manner. 

E[a] = [F'^ll' {Ql[Ro] - [Ra]i m {di m Q L [Ra]}} . (106) 

We will drop the superscript Z for simplicity but any conclusion drawn below will be valid for both specific cases i.e. Z G {(2, 2), (3, 1)}. The 
normalisation constant which acts also as inverse of associated Fisher matrix F LL i can be written as: 

F w = {(E L )(E L ,)) - {(E L )){(E L ,)) = h,di m Q L [Ra\di m Q h [Ra\) - ^(di m Q L [Ra]){di m Q L '[Ra]}. (107) 

The construction of F LL i is equivalent to the calculation presented for the case of R — C _1 . For one-point estimator we similarly can write F R — 
^2 LL , F£" L , . The optimal weighting can be replaced by arbitrary weighting. As a spacial case we can also use no weighting at all R = I. Although the 
estimator remains unbiased the scatter however increases as the estimator is no longer optimal. Use of arbitrary weights makes the estimator equivalent 
to a PCL estimator. 
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7.4 Joint Estimation of Multiple trispectra 

It may be of interest to estimate several trispectra jointly. In such scenarios it is indeed important to construct a joint Fisher matrix which will 

^[^EE^^w- (108) 

XY LL' 

The estimator E^[a] is generic and it could be either E^ 3 ' 1 ^ or E^ 2 ' 2 \ Here X and Y corresponds to different trispectra of type X and Y, these could 
be e.g. primordial trispectra from various inflationary scenarios. It is possible of course to do a joint estimation of primary and secondary trispectra. 
The off-diagonal blocks of the Fisher matrix will correspond to cross-talks between various types of bispectra. Indeed principal component analysis or 
Generalised eigenmode analysis can be useful in finding how many independent components of such trispectra can be estimated from the data. 
The cross terms in the Fisher matrix elements will be of following type: 



F £' = (i) E E (-i) M (-i) M 'W(i)] x [if;'iV(i')] y 

ST.S'T' (all Im.l'm') 

L h S \ ( L l d 8 \( l' a l' b S' \ ( L' l' d S' 
m a m b T J \ M m d -T J \ m' a m' b T' J \ M' m' d -T 

X 'LM,L> M'^l a m a ,l' a m'J^l b m b ,l' b m' b ^l a m c ,l' e m' c + ^^LM ,l' a m' a m a , V M' ^l b m b .ijm^ ^l c m c ,l' c m' c j ■ (109) 

The expression displayed above is valid only for E^- 3 ' 1 ', exactly similar results holds for the other estimator _B' 3,1 '. For X — Y we recover the results 
presented in previous section for independent estimates. As before we recover the usual result for one-point estimator for Q 4 from the Fisher matrix of 

(3 1) (2 2) 

Q L ' or Q L ' , with corresponding estimator modified accordingly. 

f xy = J2 F *^' = ^[F-r y ^H. (HO) 

LL' XY 

A joint estimation can provide clues to cross-contamination from different sources of trispectra. It also provide information about the level of 
degeneracy involved in such estimates. 



8 CONCLUSIONS 

In the near future the all-sky Planck satellite will complete mapping the CMB sky in unprecedented detail, covering a huge frequency range. The 
cosmological community will have the opportunity to use the resulting data to constrain available theoretical models. While the power spectrum provides 
the bulk of the information, going beyond this level will lift degeneracies among various early universe scenarios which otherwise have near identical 
power spectra. The higher-order spectra are the harmonic transforms of multi-point correlation functions, which contain information which can be difficult 
to extract using conventional techniques. This is related to their complicated response to inhomogeneous noise and sky coverage. A practical advance is 
to form collapsed two-point statistics, constructed from higher-order correlations, which can be extracted using conventional power- spectrum estimation 
methods. 

We have specifically studied and developed three diff erent types of estim ators which can be employed to analyse these power spectra associated 
with higher-order statistics . The MASTER-based approach jHivon et alj|200lh is typically employed to estimate pseudo-C;s from the masked sky in the 
presence of noise; also see Efstathiou (2004, 2006). These are unbiased estimators but the associated variances and scatter can be estimated analytically 

(2 1) 

with very few simplifying assumptions. We extend these estimators to study higher-order correlation functions. We develop estimators for C\ ' for the 

( 3 1 ) (2 2^ 

skew-spectrum (3-point) as well as C ; ' and C} ' which are power spectra of fields related to the trispectrum, or kurt-spectrum (4-point). The removal 
of the Gaussian contribution is achieved by applying the same es timators to a set of Gaussian simula ti ons with identical power s pectra and subtracti ng. 

As a next step we generalised the estimators employed by Yadav, Ko matsu.& Wandeltl d2007l) ; lYadav & Wandeltl §008); Yada v et al. I j2008h and 
others to study the kurt-spectrum. These methods are computationally expensive and can be implemented using a Monte-Carlo pipeline which can 
generate 3D maps from the cut-sky harmonics using radial integrations of a target theoretical model. The Monte-Carlo generation of 3D maps is the most 
computationally expensive part and dominates the calculation. The technique nevertheless has been used extensively, as it remains highly parallelisable 
and is optimal in the presence of homogeneous noise and near all-sky coverage. The corrective terms involves linear and quadratic terms for lack of 
spherical symmetry due to inhomogeneous noise and partial sky coverage. These terms can be computed using a Monte-Carlo chain. We also showed 
that the radial integral involved at the three-point analysis needs to extended to a double-integral for the trispectrum. The speed of this analysis depends 
on how fast we can generate non-Gaussian maps. It is possible also to use the same formalism to study cross-contamination from other sources such as 
point sources or other secondaries while also determining primordial non-Gaussianity. The analysis also allows us to compute the overlap or degeneracy 
among various theoretical models for primordial non-Gaussianity. 

Finally we presented the analysis fo r the case of estimators which a re completely optimal even in the presence of inhomogeneous noise and arbitrary 
sky coverage dSmifh & Zaldarriagall2006h . Extending previous work by iMunshi & Heavens! d2009l) which concentrated only on the skew-spectrum, we 
showed how to generalise to the trispectrum. This involves finding a fast method to construct and invert the covariance matrix Ci m i' m i in multipole 
space. In most practical circumstances it is possible only to find an approximation to the exact covariance matrix, and to cover this we present analysis 
for an approximate matrix which can be used as instead of C~ . This makes the method marginally suboptimal but it remains unbiased. The four-point 
correlation function also takes contributions which are purely Gaussian in nature. The subtraction of these contributions is again simplified by the use of 
Gaussian Monte-Carlo maps with the same power spectrum. A Fisher analysis was presented for the construction of the error covariance matrix, allowing 
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joint estimation of trispectra contributions from various sources, primaries or secondaries. Such a joint estimation give us fundamental limits on how 
many sources of non-Gaussianity can be jointly estimated from a specific experimental set up. A more detailed Karhuenen-Loeve eigenmode analysis 
will be presented elsewhere. The detection of non-primordial effects, such as weak lens ing of the CMB, the kinetic SZ effect, the Ostriker-Vishniac effect 
and the thermal SZ effect can provide valuable additional cosmological information dRiquelmel & Spe rgel 2006). The detection of CMB lensing at the 
level of the bispectrum needs external data sets which trace large-scale structures, but the trispectrum offers an internal detection, albeit with reduced S/N, 
providing a valuable consistency check. 

At the level of the bispectrum, primordial non-Gaussianity can for many models be described by a single parameter /jv l ■ The two degenerate kurt- 
spectra (power spectra related to tri-spectrum) we have studied at the 4-point level require typically two parameters such as /ml and gNL- Use of the 
two power spectra will enable us to put separate constraints on Jnl and gNL without using information from lower-order analysis of bispectrum, but 
they can all be used in combination. Clearly at even higher-order more parameters will be needed to describe various parameters (/jvl, gNL, fiNL, • • • ) 
which will all be essential in describing degenerate sets of power spectra associated with multispectra at a specific level. Note we should keep in mind 
that higher-order spectra will be progressively more dominated by noise and may not provide useful information beyond a certain point. 

The power of the estimators we have constructed largely depends on finding a techniques to simulate non-Gaussian maps with a specified bispectrum 
and trispectrum. We have discussed the possibility of using our technique t o generate all-sky CMB m aps with specified lower-order spectra, i.e. power 
spectrum, bispectra and trispectra. These will generalise previous results bv lSmith & Zaldarriad f2006>) which can guarantee a specific form at bispectrum 
level. 

While we have primarily focused on CMB studies, our estimators can also be useful in other areas e.g. studies involving 21cm studies, near all-sky 
redshift surveys and weak lensing surveys. The estimators described here can be useful for testing theories for primordial and/or gravity-induced secondary 
trispectrum using such diverse data sets. In the near future all-sky CMB experiments such as Planck can provide maps covering a huge frequency range 
and near-all sky coverage. The estimators described here can play a valuable role in analysing such maps. 
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APPENDIX A: 3J SYMBOLS 



Y (2/3 + l)( h h h ) ( Zl , h , 1 ) =6 mim ,5 m2m , (Al) 
/-^ v I mi m 2 mz j \ m\ m 2 ml 1 i 2 2 



y- ( h h h \( h h 1' 3 \ = S hl'Jm 3m ' 3 
I mi m 2 m 3 J \ mi m 2 ra 3 J 2/3 + 1 



